home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / fprompt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  7.1 KB  |  295 lines
  1. /* 
  2.    Find Prompt peaks:
  3. */
  4.  
  5. #include <stdio.h>
  6. #include <spec.h>
  7. #define TMPFILE "global.def"
  8.  
  9. float *spc,*err,*tim;
  10. int range=1024,
  11.     tri=0,
  12.     dev=50;
  13.  
  14. help()
  15. {
  16.   printf("finding prompt peaks in file (first approximation)\n");
  17.   printf("fprompt file [options]\n");
  18.   printf("print out prompt peak positions\n");
  19.   printf("options:\n");
  20.   printf("  -range  n  specifies the length of each spectrum (%d)\n",range);
  21.   printf("  -dev n     specifies the maximum deviation of the prompt (%d)\n",dev);
  22.   printf("  -tri [1/2] triangular spectrum beginning with minimum: 1\n");
  23.   printf("                                                peak:    2\n");
  24.   printf("  -get file  get old spectra definition file (for reading\n");
  25.   printf("                              the spectrum phase and parity)\n");
  26.   printf("  -def file  specifies another definition file for output\n\n");
  27.   printf("the default output file is %s\n",TMPFILE);
  28.   printf("This file is read by the RVT routine and is organized as follows:\n");
  29.   printf("Each line consists of 4 integers.\n");
  30.   printf("The first specifies the Time 0 position.\n");
  31.   printf("The second can be +1 -1 +2 or -2 and specifies, if it is\n");
  32.   printf("     to be mirrored (-) and if it is a 0 degree (1)\n");
  33.   printf("     or 90 degree (2) spectrum\n");
  34.   printf("The last two numbers are used to specify a range for\n");
  35.   printf("the RVT routine and the exponential fit.\n");
  36.   printf("\n(C) Rainer Kowallik\n");
  37.   exit(0);
  38. }
  39.  
  40. main(argc,argv)
  41. int argc;
  42. char *argv[];
  43. {
  44. int xx,yy,n,m,i,max,peakptr,peaks[30],phase[30],xend[30];
  45. char z[80],defnam[80];
  46. float x,y;
  47. FILE *fp;
  48.  
  49.    spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  50.    err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  51.    tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  52.  
  53.    strcpy(defnam,TMPFILE);
  54.  
  55.    guessphase(phase,tri);
  56.    getphase(defnam,phase,&tri,&range);
  57.  
  58.    if(checkopt(argc,argv,"-range",z)) {
  59.       range=atoi(z);
  60.       if(tri>0) range = 2 * range;
  61.    }
  62.    if(checkopt(argc,argv,"-dev",z)) dev=atoi(z);
  63.    if(checkopt(argc,argv,"-tri",z)) {
  64.       tri=atoi(z);
  65.       if(tri>2) {
  66.          printf("don't fool me , or I will fool you too !\n");
  67.          exit(0);
  68.       }
  69.       if(tri>0) range = 2 * range;
  70.       guessphase(phase,tri);
  71.    }
  72.    if(checkopt(argc,argv,"-def",z)) strcpy(defnam,z);
  73.  
  74.    if(checkopt(argc,argv,"-get",z)) getphase(z,phase,&tri,&range);
  75.  
  76.    max=readspec(argv[1],spc,err,tim,z);
  77.  
  78.    /* finding FIRST peak */
  79.    
  80.    yy=0; xx=0;
  81.    for(n=0;n<range;n++) {
  82.       if(spc[n]>yy) {
  83.          yy=spc[n];
  84.          xx=n;
  85.       }
  86.    }
  87.    peakptr=0;
  88.    peaks[peakptr++]=xx;
  89.    if(tri>0) peaks[peakptr++]=xx;
  90.  
  91.    /* finding all other peaks */
  92.  
  93.    n=xx+range-dev;
  94.    while(n<max) {
  95.       xx=n; yy=0;
  96.       for(i=0;i<(2*dev);i++) {
  97.          if(spc[n]>yy) {
  98.             yy=spc[n];
  99.             xx=n;
  100.          }
  101.          n = n + 1;
  102.       }
  103.       peaks[peakptr++]=xx;
  104.       if(tri>0) peaks[peakptr++]=xx;
  105.  
  106.       n=xx+range-dev;
  107.    }
  108.  
  109.    /* --------------------------------------------------------------------
  110.      finding arithmetic mean of peak heights. Used to figure out, if the
  111.      peak is real or not
  112.    --------------------------------------------------------------------- */
  113.    
  114.    xx=0;
  115.    for(i=0;i<peakptr;i++) {
  116.       xx=xx+spc[peaks[i]];
  117.    }
  118.    xx = xx/peakptr;
  119.    xx = xx / 10;
  120.  
  121.    /* --------------------------------------------------------------------
  122.       removing random "peaks"
  123.       -------------------------------------------------------------------- */
  124.   
  125.    i=0;
  126.    while(i<peakptr) {
  127.       if(spc[peaks[i]] < xx) {
  128.      for(n=i+1;n<peakptr;n++) peaks[n-1] = peaks[n];
  129.      peakptr = peakptr - 1;
  130.       } else {
  131.      i=i+1;
  132.       }
  133.    }
  134.  
  135.    /* ---------------------------------------
  136.       finding end of usable data for spectra
  137.       --------------------------------------- */
  138.  
  139.    findusable(peaks,range,tri,phase,xend,peakptr);   
  140.    n=xend[0];
  141.    for(i=0;i<peakptr;i++) {
  142.       if(xend[i]<n) n=xend[i];
  143.    }
  144.    xend[peakptr-1] = n;                   /* minimum must be the last entry */
  145.  
  146.  
  147.    /* ---------------------------------------------------
  148.       writing data 
  149.    ------------------------------------------------------ */
  150.  
  151.    fp=fopen(defnam,"w");
  152.    for(i=0;i<peakptr;i++) {
  153.       fprintf(fp,"%d   %2d   %d   %d\n",peaks[i],phase[i],5,xend[i]);
  154.    }
  155.    fclose(fp);
  156.    printf("\ndata written to %s\n now start editor on this file !\n",TMPFILE);
  157.    free(spc); free(err); free(tim);
  158.    exit(0);
  159. }
  160.  
  161. iabs(x)
  162. int x;
  163. {
  164. if(x>=0) return(x);
  165. return(-x);
  166. }
  167.  
  168.  
  169. /* -----------------------------------------------
  170.    reading spectra definition file
  171.    ----------------------------------------------- */
  172. getphase(name,phase,tri,range)
  173. char name[];
  174. int phase[], *tri, *range;
  175. {
  176. int i,n,m,p[30];
  177. FILE *fp;
  178.  
  179.    fp=fopen(name,"r");
  180.    if(fp==NULL) {
  181.       return(0);
  182.    }
  183.    i=0; *tri=0;
  184.    while(!feof(fp)) {
  185.       fscanf(fp,"%d %d %d %d\n",&p[i],&phase[i],&n,&n);
  186.       i=i+1;
  187.    }
  188.    fclose(fp);
  189.    for(n=0;n<i;n++) if(phase[i]<0) *tri=2;   /* triangular spectrum ? */
  190.    if(phase[0]<0) *tri=1;        /* triangular, beginning with minimum */
  191.    n=p[0];
  192.    if(p[1]>(n+5)) {
  193.       n = iabs(n - p[1]);
  194.    } else {
  195.       n = iabs(n - p[2]);
  196.    }
  197.    i=1;
  198.    while(i<n) i= i << 1;
  199.    m = i >> 1;
  200.    if(iabs(n-m) < iabs(n-i)) {
  201.       n = m;
  202.    } else {
  203.       n = i;
  204.    }
  205.    *range = n;
  206. }
  207.  
  208. guessphase(phase,tri)
  209. int phase[], tri;
  210. {
  211. int i,n,m;
  212.  
  213.    n=1;
  214.    for(i=0;i<30;i++) {
  215.       n= -1 * n;
  216.       switch(tri) {
  217.       case 0:
  218.          if(n>0) phase[i]=1;
  219.          if(n<0) phase[i]=2;
  220.          break;
  221.       case 1:
  222.          if(n>0) {
  223.             phase[i++] = -1;
  224.             phase[i]=1;
  225.          }
  226.          if(n<0) {
  227.             phase[i++]= -2;
  228.             phase[i]=2;
  229.          }
  230.          break;
  231.       case 2:
  232.          if(n>0) {
  233.             phase[i++]=1;
  234.             phase[i]= -1;
  235.          }
  236.          if(n<0) {
  237.             phase[i++]=2;
  238.             phase[i]= -2;
  239.          }
  240.          break;
  241.       }
  242.    }
  243. }
  244.  
  245. /* ------------------------------------------------
  246.    finding usable data area
  247.    ------------------------------------------------ */
  248. findusable(peaks,range,tri,phase,xend,peakptr)
  249. int peaks[],phase[],xend[],
  250.     range,tri,peakptr;
  251. {
  252. int i,n,m,start,end;
  253.  
  254.    start = range / 2;
  255.    end = range;
  256.    if(tri>0) {
  257.       start = range / 4;
  258.       end = range / 2;
  259.    }
  260.    end = end - 3;
  261.    for(n=0;n<peakptr;n++) {  /* scanning spectra */
  262.       xend[n] = end;
  263.       if(phase[n] > 0) {    /* positive time axis */
  264.          m = peaks[n] + start;
  265.          for(i=start;i<end;i++) {
  266.             m = m + 1;
  267.             if(!outrange(spc[m],err[m],spc[m+1],err[m+1])) continue;
  268.             if(!outrange(spc[m],err[m],spc[m+2],err[m+2])) continue;
  269.             if(!outrange(spc[m-1],err[m-1],spc[m+1],err[m+1])) continue;
  270.             xend[n] = i - 2 ;
  271.          }
  272.       } else {
  273.          m = peaks[n] - start;
  274.          for(i=start;i<end;i++) {
  275.             m = m - 1;
  276.             if(!outrange(spc[m],err[m],spc[m-1],err[m-1])) continue;
  277.             if(!outrange(spc[m],err[m],spc[m-2],err[m-2])) continue;
  278.             if(!outrange(spc[m+1],err[m+1],spc[m-1],err[m-1])) continue;
  279.             xend[n] = i - 2 ;
  280.          }
  281.       }
  282.    }
  283. }
  284.  
  285. outrange(x1,e1,x2,e2)
  286. float x1,e1,x2,e2;
  287. {
  288. float a,b,c;
  289.  
  290.    a = (float) fabs((double)(x1-x2));
  291.    b = 3 * (e1 + e2);
  292.    if(a>b) return(TRUE);
  293.    return(FALSE);
  294. }
  295.